
from math import radians, cos, sin, asin, sqrt, pow
'''
# Douglas - Peuker算法(DP算法)
# 过程如下:
# 1、连接曲线首尾两点A、B；
# 2、依次计算曲线上所有点到A、B两点所在曲线的（直线）距离；
# 3、计算最大距离D，如果D小于阈值threshold, 则去掉曲线上出A、B外的所有点；如果D大于阈值threshold, 则把曲线以最大距离分割成两段；
# 4、对所有曲线分段重复1 - 3
# 步骤，知道所有D均小于阈值。即完成抽稀。

使用的抽稀算法，注意这里计算的距离是利用经纬度与地球信息计算的真实距离
'''




def verticalDist(a, b, c):
    '''

    计算a点 到 b 与 c 点组成直线的距离

    :param a:
    :param b:
    :param c: 目标点
    :return:
    '''

    S, L = _triangleArea(a, b, c)
    return 2 * S / L


def _triangleArea(pointA, pointB, pointC):
    '''
      三角形面积。海伦-秦九韶公式
      s = 开根(p(p-a)(p-b)(p-c))
    '''
    a = _getDistance(pointA, pointB)
    b = _getDistance(pointA, pointC)
    c = _getDistance(pointB, pointC)
    p = (a + b + c) / 2
    return sqrt(p*(p-a)*(p-b)*(p-c)), c


def _getDistance(pointA, pointB):
    '''
      根据两点经纬度，获取两点间的距离
      point = [lon, lat]
    '''
    lng1, lat1, lng2, lat2 = map(
        radians, [pointA[0], pointA[1], pointB[0], pointB[1]])  # 经纬度转换成弧度
    dlon = lng2 - lng1
    dlat = lat2 - lat1
    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
    distance = 2 * asin(sqrt(a)) * 6371 * 1000  # 地球平均半径，6371km
    distance = round(distance, 3)
    return distance


def point2LineDistance(point_a, point_b, point_c):
    """
    计算点a到点b c所在直线的距离， 并非地球的球面距离。

    这里只能从理解的角度将这个距离与地球经纬度距离相关联！

    point, point_list[0], point_list[-1]
    :param point_a: 中间点 [ ]
    :param point_b: 线段起点 []
    :param point_c: 线段终点 []
    :return:
    """
    # 首先计算b c 所在直线的斜率和截距
    if point_b[0] == point_c[0]:
        return 9999999
    slope = (point_b[1] - point_c[1]) / (point_b[0] - point_c[0])
    intercept = point_b[1] - slope * point_b[0]

    # 计算点a到b c所在直线的距离
    distance = abs( slope * point_a[0] - point_a[1] + intercept) / sqrt(1 + pow(slope, 2) )
    return distance


class DouglasPeuker(object):
    def __init__(self , THRESHOLD):
        self.threshold = THRESHOLD
        self.qualify_list = list() # 最终返回的抽样后的点
        self.disqualify_list = list() # 计算完成时，这个列表为空

    def diluting(self, point_list):
        """
        抽稀过程
        :param point_list:[[], [], [], ....]二维点列表
        :return:
        """
        if len(point_list) < 3:
            self.qualify_list.extend(point_list[::-1])
        else:
            # 找到与收尾两点连线距离最大的点
            max_distance_index, max_distance = 0, 0
            for index, point in enumerate(point_list):
                if index in [0, len(point_list) - 1]: # 定义边界
                    continue

                # distance = point2LineDistance(point, point_list[0], point_list[-1]) # 计算直线距离
                distance = verticalDist(point, point_list[0], point_list[-1]) # 单个点到线段首尾的距离
                # 计算首尾两端 最大距离的点
                if distance > max_distance:
                    max_distance_index = index
                    max_distance = distance

            # 若最大距离小于阈值，则去掉所有中间点。
            if max_distance < self.threshold:
                self.qualify_list.append(point_list[-1])
                self.qualify_list.append(point_list[0])
            else:
                # 将曲线按最大距离的点分割成两段
                sequence_a = point_list[:max_distance_index]
                sequence_b = point_list[max_distance_index:]

                for sequence in [sequence_a, sequence_b]:
                    if len(sequence) < 3 and sequence == sequence_b:
                        self.qualify_list.extend(sequence[::-1])
                    else:
                        self.disqualify_list.append(sequence)

    def main(self, point_list):
        '''
        :param point_list:  [[long , lat] , [...] , ...]  经度维度list
        :return:
        '''

        self.diluting(point_list) # 首先对这条线进行一次抽稀计算，若 len(self.disqualify_list == 0 ，则一次计算完成，否则需要继续抽稀
        while len(self.disqualify_list) > 0:
            self.diluting(self.disqualify_list.pop()) # 这里的被pop()抛出的对象会被self.diluting() 接收
        return self.qualify_list